Generalized harmonic spatial coordinates and hyperbolic shift conditions 
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We propose a generalization of the condition for harmonic spatial coordinates analogous to the 
generalization of the harmonic time slices introduced by Bona et at, and closely related to dynamic 
shift conditions recently proposed by Lindblom and Scheel, and Bona and Palenzuela. These gen- 
eralized harmonic spatial coordinates imply a condition for the shift vector that has the form of an 
evolution equation for the shift components. We find that in order to decouple the slicing condition 
from the evolution equation for the shift it is necessary to use a rescaled shift vector. The initial 
form of the generalized harmonic shift condition is not spatially covariant, but we propose a simple 
way to make it fully covariant so that it can be used in coordinate systems other than Cartesian. 
We also analyze the effect of the shift condition proposed here on the hyperbolicity of the evolution 
equations of general relativity in 1+1 dimensions and 3+1 spherical symmetry, and study the pos- 
sible development of blow-ups. Finally, we perform a series of numerical experiments to illustrate 
the behavior of this shift condition. 
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I. INTRODUCTION 

When one studies the time evolution of the gravita- 
tional field in general relativity, a good choice of coor- 
dinates (a "gauge" choice) can make the difference be- 
tween finding a well behaved solution for a large portion 
of the spacetime, or running into a coordinate (or physi- 
cal) singularity in a finite coordinate time, which would 
not allow a numerical evolution to continue any further. 
In the 3+1 formulation, the choice of the time coordi- 
nate is related with the lapse function, while the choice 
of the spatial coordinates is related to the shift vector. 
Many different ways to choose the lapse and the shift 
have been proposed and used in numerical simulations in 
the past (see for example the pioneering papers of Smarr 
and York 

TO) 

Some gauge choices involve solving ellip- 
tic equations, while others involve solving evolution type 
equations, which may or may not be hyperbolic in char- 
acter. Recently, hyperbolic coordinate conditions have 
become a focus of attention, as they in principle allow 
one to write the full set of dynamical equations as a 
well-posed system 0, 0, IE El 0- while at the same 
time being both easier to implement and considerably 
less computationally expensive than elliptic conditions. 

The classic example of hyperbolic coordinate condi- 
tions are the so-called harmonic coordinates, which are 
defined by asking for the wave operator acting on the 
coordinate functions to vanish. Harmonic coordinate 
conditions have the important property of allowing the 
Einstein field equations to be written as a series of wave 
equations (with non-linear source terms) for the metric 
coefficients g^ v . Because of this, these conditions were 
used to prove the first theorems on the existence of so- 
lutions to the Einstein equations 9]. This property of 
transforming the Einstein equations into wave equations 
could in principle also be seen as an important advantage 



in the numerical integration of these equations. Still, 
with few exceptions (see for example 0, 0, H^). full 
harmonic coordinates have traditionally not been used 
in numerical relativity, though harmonic time slices have 
been advocated and used in some cases [lj, Hj, Ha lilil • 
The reason for this is two-fold: In the first place, har- 
monic coordinates are rather restrictive, and formula- 
tions of the Einstein equations for numerical relativity 
are usually written in a way that allows the gauge free- 
dom to remain explicit so it can be used to control certain 
aspects of the evolution (avoid singularities, enforce sym- 
metries, reduce shear, etc.). Also, in the particular case 
of a harmonic time coordinate, it has been shown that 
the space-like foliation avoids focusing singularities only 
marginally, and is therefore not a good choice in many 
cases |g, [13J, llll llSj- Of course, it can be argued that 
any coordinate choice is harmonic if one does not ask for 
the wave operator acting on the coordinate functions to 
be zero, but instead to be equal to a known function of 
spacetime (a "gauge source function"). This is certainly 
true, but of little use in real life numerical simulations 
where there is no way to know a priori what is a conve- 
nient choice for these gauge source functions (but see 01 
for some suggestions that seem to work well in practice). 

Nevertheless, the fact that the use of harmonic co- 
ordinates allows the field equations to be written in 
strongly hyperbolic form makes one immediately ask if 
there might be simple generalizations of the harmonic 
conditions that will still allow the field equations to be 
written in strongly hyperbolic form, while at the same 
time retaining a useful degree of gauge freedom. That 
this is indeed the case was first shown for the particular 
case of a harmonic time coordinate by Bona et al. in |19| , 
where a strongly hyperbolic reformulation of the Einstein 
evolution equations was constructed using a generalized 
harmonic slicing condition which is usually referred to as 
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the Bona-Masso slicing condition. It includes as particu- 
lar cases several choices that had been used in numerical 
simulations from the early 90's wit h g ood results, such 
as for example the "1+log" slicing [2(J,E]]- In fact, the 
Bona-Masso slicing condition was motivated precisely to 
include such empirically tested conditions in a strongly 
hyperbolic formulation of the Einstein equations. 

In this paper we want to follow a similar approach and 
propose a generalization of the harmonic spatial coordi- 
nate condition. We will show how this allows us to obtain 
a hyperbolic shift condition that is very closely related 
to conditions already proposed in the literature, most 
notably the shift conditions recently introduced byLind- 
blom and Scheel , and by Bona and Palenzuela |22j | (in 
fact, under some specific circumstances, one finds that 
the shift condition proposed here becomes a particular 
case of those of Refs. p and [22|). 

This paper is organized as follows. In Sec. [pj we dis- 
cuss the standard harmonic coordinates and write them 
as evolution equations for the lapse and shift. We also in- 
troduce a rescaled shift vector that allows one to decouple 
the lapse and shift equations. Section ITTT1 generalizes the 
condition for spatial harmonic coordinates, and Sec. IIVI 
discusses the interpretation of this condition in curvilin- 
ear coordinate systems. In Sec. we describe the con- 
cept of hyperbolicity and the source criteria for avoiding 
blow-ups. Section IVII studies the generalized harmonic 
shift condition in the case of 1+1 dimensions, analyzing 
its hyperbolicity properties, the possible appearance of 
blow-ups (gauge shocks), and also the behavior of this 
shift condition in numerical simulations. In Sec. IVIll we 
repeat the same type of analysis for spherical symmetry 
and also present results from numerical simulations. We 
conclude in Sec. I Villi Finally, Appendix El shows a for- 
mal derivation of the generalized harmonic lapse and shift 
conditions, and Appendix [51 gives general expressions for 
the 4-Christoffcl symbols in terms of 3+1 quantities. 

II. HARMONIC COORDINATES 

Let us consider four scalar coordinate functions <j) a de- 
fined on a given background spacetime. The condition 
for these coordinates to be harmonic is simply 

Utj> a :=5""V M V^ Q =0, (2.1) 

with the spacetime metric tensor. 

Let us further assume that (f>° is such that its level 
surfaces are space-like. In that case, 0° can be identified 
with a global time function. If we define the lapse func- 
tion a as the interval of proper time when going from the 
hypersurface cjp = t to the hypersurface (f>° = t + dt along 
the normal direction, then it is easy to show that a will 
be given in terms of (jP as 

a= (-V0°-V0°)" 1/2 . (2.2) 

The definition of the shift vector is somewhat more 
involved. We start by defining three scalar functions (3 a 



such that when we move from a given level surface of (jP 
to the next following the normal direction, the change in 
the spatial coordinate functions <f> a is given by 

$ +dt = ft- (3 a d<tP , (2.3) 
from which one can easily find 

(3 a = - a (n- V0 a ) , (2.4) 
with n the unit normal vector to the hypersurface <f)° = t, 

n = -aV(j) 1 (2.5) 

and where the minus sign is there to guarantee that n 
is future pointing. Thus defined, the (3 a are scalars, but 
we can use them to define a vector /3 by asking for its 
components in the coordinate system {<p a } to be given 
by (0,/3 a ). The vector constructed in this way is clearly 
orthogonal to n. In an arbitrary coordinate system {x^}, 
the shift components will then be given by 

0» = -a(fi.Vr)^. (2.6) 

Notice that with this definition, the shift vector is pro- 
portional to the lapse function, so that a simple rescaling 
of (fP changes the shift. This suggests that it is perhaps 
more natural to define a rescaled shift vector a in the 
following way 

„» : =2- = _ (S . vr) _. ,2.7, 

We will see below that this rescaled shift vector will be 
important when expressing the harmonic condition in 
3+1 language. 

The harmonic coordinate conditions can be simplified 
by expanding them in the coordinate system {x a = cf> a }, 
in which case they reduce to 

T" := <rr^ = , (2.8) 

where r" are the Christoffel symbols associated with the 
4-metric g^. If we now relate the coordinates {x a = <j> a } 
to the standard 3+1 coordinates, then these four equa- 
tions can be shown to become (see Appendix E] and [BJ) 

d t a = !3 a d a a-a 2 K, (2.9) 
fit/3* = (3 a d a (3 l - ad 1 a + a 2 

+— (dta - P a d a a + a 2 K) . (2.10) 
a 

Here K is the trace of the extrinsic curvature, and ( 3 'r* is 
defined in terms of the three-dimensional Christoffel sym- 
bols ^ s 'Tj k , the spatial metric jij and its determinant 

7 := det 7ij by ( 3 )P := ^T) k = -d, (^7 7^) ly/l- 
Notice that in equation l|2.1(J|l we have an explicit depen- 
dency on the time derivative of the lapse function. This 
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dependency is usually not written down, as the whole last 
term of the second equation vanishes if the first equa- 
tion is assumed to hold, but we prefer to leave the de- 
pendency explicit (see for example 111! | 23lk incidentally, 
equation (|2.10() fixes a sign error in [23J, and includes a 
term missing in jllj). 

The fact that the evolution equation for the shift de- 
pends on the time derivative of the lapse is inconvenient 
if one wants to use harmonic spatial coordinates with a 
different slicing condition, say maximal slicing. It is also 
an indication that the shift itself might not be the most 
convenient function to evolve. Remarkably, it turns out 
that if we rewrite the evolution equation for the shift in 
terms of the rescaled shift a 1 — p l /a introduced above, 
then the spatial harmonic condition decouples completely 
from the evolution of the lapse. We find 

d t <7 1 =aa a d a (i l -d l a + a(a l K + < 3 'r 4 ) . (2.11) 

Therefore, if one works with a 1 instead of one can 
use harmonic spatial coordinates with an arbitrary slicing 
condition in a straightforward way. 

A final comment about equations i|2.9|) and i|2.10[l is 
in order. Equation (|2.9|) is clearly a scalar equation as 
seen in the spatial hypersurfaces. Equation l|2.10|l . on 
the other hand, is not 3-covariant, i.e. starting from 
exactly the same 3-geometry but in different coordinates, 
it will produce a different evolution for the shift vector. 
This might seem surprising since this equation is just 
the 3+1 version of the condition for spatial harmonic 
coordinates which is 4-covariant. However, there is no 
real contradiction, since changing the coordinates on the 
spatial hypersurfaces means changing the scalar functions 
(f> 1 themselves, so it should not be surprising that we get a 
different shift. We will come back to this point in Sec. llVI 
where we will propose a way to make the shift evolution 
equation fully 3-covariant. 

III. GENERALIZED HARMONIC 
COORDINATES 

In 0| i Bona et al. generalize the harmonic slicing 
condition (|2.9|l in the following way 

d t a- (3 a d a a = -a 2 f{a)K , (3.1) 

with f(a) a positive but otherwise arbitrary function of 
the lapse. This slicing condition was originally motivated 
by the Bona-Masso hyp erbolic reformulation of the Ein- 
stein equations [Til Il5l fl8l \l$L |24| , but it can in fact be 
used with any form of the 3+1 evolution equations. As 
discussed in the Bona-Masso slicing condition above 
can be shown to avoid both focusing singularities 
and gauge shocks [25[ for particular choices of /. Refer- 
ence p| also shows that condition l|3.1|l can be written in 
4-covariant form in terms of a global time function </>° as 

(^-a/nVjV^^O, (3.2) 



with a i :— 1/ f(a) — 1 and the unit normal vector 
to the spatial hypersurfaces defined in (|2.5|) . Here we 
will introduce an analogous generalization of the spatial 
harmonic coordinates {</>'}. That is, we propose the fol- 
lowing spatial gauge condition 

-a h n^n v )V^^ 1 = , (3.3) 

where is still the unit normal to the spatial hypersur- 
faces, but now at ■= 1/h — 1, with h(a, (3" 1 ) a scalar func- 
tion that can in principle depend on both the lapse and 
shift (we will see below that the shift dependence is in fact 
not convenient). In the coordinate system {x M = (f)^}, 
condition (|3.3|l becomes 

(f-a k nV)r|,„=0. (3.4) 

Expressing the 4-metric and normal vector in terms of 
3+1 variables, the last equation becomes 

r< - 2/3™rJ„ + (3 m P n r l mn = a 2 h 7 mn r l mn . (3.5) 

Notice that on the right hand side of this equation ap- 
pears the contraction r ) mn T l mn which should not be con- 
fused with r ; :— g^ LV T l ^ v . Inserting now the expressions 
for the T l mn in terms of 3+1 quantities we obtain (see 
also Appendix |B|) 

Rl 

d t p l = (3 m d m (3 l - ad l a + — (d t a - f3 m d m a) 

a 

+ a 2 h(^K + WT l \ . (3.6) 

This is to be compared with equation (|2.10(1 of the previ- 
ous section. Notice that again we find that the evolution 
equation for the shift is coupled to that of the lapse. In 
the same way as before, we can decouple the shift evolu- 
tion equation by writing it in terms of the rescaled shift 
a 1 = (3 l /a. We find 

t a l =aa m d m a l - d l a + ah(a l K + ^T l ^j , (3.7) 

which is to be compared with (|2.11|) . This is the final 
form of the condition for generalized harmonic spatial 
coordinates, and we will refer to this condition simply as 
the "generalized harmonic shift" (but see Sec. II VI below 
where the condition is somewhat modified to make it fully 
3-covariant). 

At this point it is important to discuss the relation that 
the shift condition l|3.6|l has with the conditions recently 
proposed by Lindblom and Scheel @, and by Bona and 
Palenzuela |22|. It is not difficult to see that by choos- 
ing the free parameters in these references appropriately, 
one can in fact recover condition (|3.6[) . but only provided 
one also takes the lapse to evolve via the Bona-Masso 
slicing condition (|3.1|l and takes / = h. If, on the other 
hand, one uses a different slicing condition (say maximal 
slicing), or uses the Bona-Masso slicing condition with 
/ 7^ h, then this is no longer the case and the shift con- 
dition proposed here will differ from those of Refs. § 
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and [22| • This is a crucial point, and shows the impor- 
tance of rescaling the shift in order to decouple its evo- 
lution equation from the time derivative of the lapse. 

In the following sections we will study this shift condi- 
tion. Wc will first discuss the issue of the interpretation 
of the generalized harmonic shift condition for curvilinear 
coordinates in Sec. II VI Later, in Sec.|V]we will introduce 
the concept of hyperbolicity, and a criteria for avoiding 
blow-ups in the solutions of strongly hyperbolic systems 
of equations. Finally, in Sections IVII and IVIII we will 
consider the special cases of 1+1 dimensional relativity 
and spherical symmetry. In each case we will analyze 
the hyperbolicity properties of the full system of equa- 
tions including the generalized harmonic shift condition, 
study the possible development of blow-ups, and present 
a series of numerical examples. 



IV. CURVILINEAR VERSUS CARTESIAN 
COORDINATES 

We have already mentioned that the harmonic shift 
condition (|2.11|) . and its generalization (|3.7|) . are in fact 
not covariant with respect to changes in the spatial co- 
ordinates. That is, starting from exactly the same 3- 
geometry but with different spatial coordinates we will 
get a different evolution of the shift vector. In particular, 
for curvilinear systems of coordinates one could find that 
even starting from a flat slice of Minkowski spacetime we 
would still have non-trivial shift evolution driven by the 
fact that the initial ( 3 )r l do not vanish (i.e. the spatial 
curvilinear coordinates are not 3-harmonic). Worse still, 
in many cases it can happen that the ( 3 )r* of flat space 
are not only non-zero but are also singular, as is the case 
with spherical coordinates for which ( 3 )r r is of order 1/r. 
One may also find that in physical systems that have a 
specific symmetry the shift evolution will break the sym- 
metry because of the properties of some of the An 
example of this are again spherical coordinates for which 
one finds that ( 3 )r e ^ 0, so a e will evolve away from zero 
even for a spherically symmetric system. 

The question then arises how to interpret the harmonic 
shift condition in a general coordinate system, and in 
particular how to make sure that we do not run into 
pathological situations like those described above. Our 
proposal for resolving this issue is to always apply the 
generalized harmonic shift condition in a coordinate sys- 
tem that is topologically Cartesian. Of course, if one 
has a situation that has a specific symmetry, one would 
like to work with a coordinate system that is adapted to 
that symmetry. We therefore need to transform condi- 
tion (|3.7|l from Cartesian coordinates to our curvilinear 
coordinates, but taking into account the fact that the 
condition is not covariant. 

Let us denote by {x a } our reference topologically 
Cartesian coordinates, and by {x 1 } the general curvi- 
linear coordinates. If we assume that condition l|3.7|) is 



satisfied for the original coordinates {x a } we will have 
d t a a =aa l d- b a a -d a a + ah(a a K + ^T a ) . (4.1) 



In order to transform this expression we will use the fact 
that with respect to the 3-geometry a % behaves like a 
vector, while a and K behave as scalars. Remembering 
now that the Christoffel symbols transform as 



^v) k = (d- a x l djx 1 d k x E ) ( 3 'rf g + 



jk > 



(4.2) 



with Fj k := daX 1 djdkX a , we find that in the curvilinear 
coordinate system equation 1)4. If) becomes 



dto = aa m d m a 



aa a i< — a a 



+ ah(a l K + ^T l -7""^™) • (4.3) 

By rearranging some terms, the shift condition can finally 
be written in the more convenient form 



d t o l = 



l V m cr' - V'a + ah a l K 



+ a (ft/7"" 1 - (J m a n ) A l r 



(4.4) 



with A l mn := ( 3 'T l mn — Ff nn . The last expression is in 
fact 3-covariant, as one can readily verify that the A l mn 
transform as the components of a 3-tensor. But the price 
we have paid is that we have chosen a privileged coordi- 
nate system to be used as a reference in order to define 
F l mn . It is clear that for the original coordinates {x a } 
the condition above reduces to what we had before since 
F l mn vanishes. We will consider the case of spherical co- 
ordinates in Sec. IVIII below. 

In practice, one can use the fact that for flat space 
in Cartesian coordinates the Christoffel symbols vanish, 
which implies 



F 



(3) 



Hat 



so that 



(3) pi _ (3) pi 



Hal 



(4.5) 



(4.6) 



V. HYPERBOLICITY AND SHOCKS 
A. Hyperbolic systems 

The concept of hyperbolicity is of fundamental im- 
portance in the study of the evolution equations asso- 
ciated with a Cauchy problem as the initial value prob- 
lem for strongly or symmetric hyperbolic systems can 
be shown to be well-posed (though the well-posedness 
of strongly hyperbolic systems requires that some addi- 
tional smoothness conditions are verified) . In the follow- 
ing we will concentrate on one-dimensional systems, for 
which the distinction between strongly and symmetric 
hyperbolic systems does not arise. 
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Following 26] , we will consider quasi-linear systems of 
evolution equations that can be split into two subsystems 
of the form 

d t u = M(u) v , (5.1) 
d t v + A(u) d x v = q v (u, v) . (5.2) 

Here it and v are n and m dimensional vector-valued 
functions, and M and A are n x n and m x m matrices, 
respectively. In addition we demand that the w's are re- 
lated to either time or space derivatives of the u's. This 
implies that derivatives of the it's can always be substi- 
tuted for w's and hence treated as source terms. 

The system of equations above will be hyperbolic if 
the matrix A has m real eigenvalues A^. Furthermore, 
it will be strongly hyperbolic if it has a complete set 
of eigenvectors If we denote the matrix of column 
eigenvectors by R = (£i • • - £m)) then the matrix A can 
be diagonalized as 

R" l AR = diag [Ai, • • • , X m ] = A . (5.3) 

For a strongly hyperbolic system we then define the 
eigenfields as 

w = R- 1 v. (5.4) 

By analyzing the time evolution of the eigenfields, one 
can identify mechanisms that lead to blow-ups in the so- 
lution, which in [27| have been referred to as " geom etric 
blow-up" (leading to "gradient catastrophes" |23) and 
the "ODE-mechanism" (causing "blow-ups within finite 
time"). In |26fl some of us presented blow-up avoiding 
conditions for both these mechanisms, which we called 
"indirect linear degeneracy" [2!| and the "source crite- 
ria" . In that reference it was also shown, using numerical 
examples, that the source criteria for avoiding blow-ups 
is generally the more important of the two conditions. 
Because of this, and also because of the fact that the 
true relevance of indirect linear degeneracy is not as yet 
completely clear, in this paper we will concentrate only 
on the source criteria. 



B. Source criteria for avoiding blow-ups 

An evolution variable can become infinite at a given 
point by a process of "self-increase" in the causal past of 
this point. A criteria to avoid such blow-ups for systems 
of partial differential equations of the form (|5.1|l - (|5.2(1 
was proposed by some of us in |26j | : When diagonalizing 
the evolution system for the v'b, making use of (|5.3|l and 
lfO|) . one finds 

d t w + A d x w = q w , (5.5) 

where 

q w := R" 1 ^ + [dtR- 1 + Ac^R" 1 ] v . (5.6) 



This yields an evolution system where on the left- 
hand side of i|5.5[l the different eigenfields Wi are decou- 
pled. However, in general the equations are still coupled 
through the source terms q Wi . In particular, if the origi- 
nal sources were quadratic in the v's, one obtains 

d rn 

—jj- = d t Wi + Kd x Wi = CijkWjWk + 0(w) , (5.7) 
j,k=i 

where d/dt := dt + \id x denotes the derivative along the 
corresponding characteristic. As pointed out in [26(, the 
cmwf component of the source term can be expected to 
dominate and to cause blow-ups in the solution within a 
finite time. In order to avoid these blow-ups we therefore 
demand that the coefficients cm should vanish, and we 
refer to this condition as the "source criteria" . 

It is not difficult to convince one-self that the source 
criteria is in fact not a sufficient condition for avoiding 
blow-ups, as already discussed in [2(|. However, one can 
still expect the source criteria to be a necessary condi- 
tion for avoiding blow-ups at least for small perturba- 
tions propagating with different eigenspeeds, as mixed 
terms will be suppressed when pulses moving at different 
speeds separate from each other, while the effect of the 
term cmw^ will remain as each pulse moves. If, however, 
some eigenfields Wi and Wj travel with identical or similar 
eigenspeeds, then one should also expect important con- 
tributions coming from the mixed terms WiWj. We will 
show an example later on where eliminating such mixed 
terms (in addition to the quadratic terms) indeed leads 
to further improvements. 

VI. EINSTEIN EQUATIONS IN 1+1 
DIMENSIONS 

We first consider standard general relativity in one spa- 
tial dimension (and in vacuum). Since in this paper we 
are interested precisely in studying a new shift condition, 
1+1 dimensional relativity is an ideal testing ground for 
the "gauge dynamics" which one can expect in the higher 
dimensional case. 

In the following sections we will introduce the evolu- 
tion equations and gauge conditions, and consider the 
possible formation of blow-ups associated with our gauge 
conditions. We will also present numerical simulations 
that show how the generalized harmonic shift condition 
behaves in practice. 

A. Evolution equations 

We will start from the "standard" Arnowitt-Deser- 
Misner (ADM) equations for one spatial dimension [2j| 
as formulated in |23j . In this case the u quantities consist 
of the lapse function a, the rescaled shift a :— a x , and 
the spatial metric function g := g xx . The v quantities, 
on the other hand, are given by the spatial derivatives 
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of the it's and, in addition, the unique component of the 
extrinsic curvature. That is, 

u = (a, a, g) , v = (^D a ,d a ,D g ,K^ , (6.1) 
with K := ^fgtvK = ^fg K£, and where we have defined 

D a := d x In a , L> 9 := d x In 5 , d CT := d x a . (6.2) 

Notice first that we use a rescaled extrinsic curvature, as 
this makes the evolution equations considerably simpler. 
Also, we use logarithmic spatial derivatives of a and g, 
but only the ordinary spatial derivative of the rescaled 
shift a, as the shift is allowed to change sign. 

For the evolution of the gauge variables we will use the 
Bona-Masso slicing condition (|3.1|l and the generalized 
harmonic shift condition (|3.7|1 . The equations for the u's 
are then 



d t ct 
d t a 
dtg 



aD a - 



ad„ 



ag 



aD„ 




(6.3) 
(6.4) 
(6.5) 



where / = f(a) and h = h(a,a). The evolution equa- 
tions for {D a ,d a ,D g } can be obtained directly from 
the above equations, while the evolution equation for K 
comes from the ADM equations and takes the following 
simple form 



d t K = d x 



a(aK - D a /y/g 



(6.6) 



The evolution equations for the u's can then be writ- 
ten in full conservative form dtv + d x (A v) = 0, with the 
characteristic matrix A given by 



A = 



— aa <xf/-\fg \ 

a/g —aa —ahjlg —aha/y/g I 

— 2aa —2a —aa 2a/^/g I 

a /\/9 —aa ) 



(6.7) 



This matrix has the following eigenvalues 
\{ = a(±^jTg-a) 



A'i 



(±^hf 9 -a) , 



(6.8) 
(6.9) 



with corresponding eigenfunctions (the normalization is 
chosen for convenience) 



K±D a /^f , 



(6.10) 



(l±a^/gt)K + ^gd a TVhD g /2, (6.11) 



da 

K 



which can be easily inverted to find 

— [w J + - w J _j , 

— !— (w{ + w f _ + w'l + w h _ 
2\/5 v 

-L( w h _-w h + )-a^g-(w{ 



(6.12) 
(6.13) 

(6.14) 
(6.15) 



The system is therefore strongly hyperbolic as long as 
/ > and h > 0, with the lapse and shift eigenficlds 
u>£ and w± propagating with the corresponding gauge 
speeds A£ and Aj_. 



B. Gauge shock analysis 

By analyzing quadratic source terms in the evolution 
equations of the eigenfields u>i , we now want to study the 
possible formation of blow-ups for the system of evolution 
equations of the previous section. For the lapse and shift 
eigenfields we find 



dw± 
~~dT 



dt 



iff /2 , ffh f h 

c±±±w± + c± ±± w±w± 



hhh „„fi2 



**±±± w ± 

O (w h ± w\,w h ± 



hhf 



f 



w+ + c±±±w±w± 

f 



(6.16) 



(6.17) 



In particular, we observe that in (|6.16|) no term propor- 
tional to w± 2 is present, and in the same way in (|6.17|) 

f 2 

there is no term proportional to w± . In order to ap- 
ply the source criteria we need to calculate those terms 
quadratic in Wi appearing in the sources of the evolution 
equation for Wi itself. It turns out that the cm coefficients 
have the form 



c f J£ ± cx (l-/-a/'/2) 



„hhh 
-±±± 



cx dh/da . 



(6.18) 
(6.19) 



According to the source criteria these coefficients have to 
vanish in order to avoid blow-ups. The conditions on the 
gauge functions /(a) and h(a, a) are then 



l-f-af/2 = 
dh/da = 



(6.20) 
(6.21) 



The condition 



times before 0, 



(a) has been studied many 
and its general solution is 



/(a) = 1 + const/a 2 



(6.22) 
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For h(a, <r), on the other hand, we obtain the condition 
that h can be an arbitrary function of a, but may not 
depend on er, that is, h = h(a). 

One now might wonder about the case where h is 
equal (or very close to) the function /. In that case 
the eigenfields w± and w± travel with the same (or sim- 
ilar) eigenspeeds, so mixed terms of the type w±w± in 
the sources can be expected to contribute to a blow-up. 
For this reason we have also calculated the cnj coeffi- 
cients associated to these terms. Notice, however, that 
in general the coefficients of such mixed terms are not 
invariant under rescalings of the eigenfields of the form 
Wi = Qi(a,a,g) Wi, so we have in fact done the calcula- 
tion assuming an arbitrary rescaling. We find 



Jfh 

-±±± 



hhf 
c ±±± °< 




(6.23) 



(6.24) 



One can readily verify that these coefficients vanish for 
/ = h = 1 + const/a 2 , independently of the rescaling of 
the eigenfields. This setting of / and h hence seems to 
be an optimal choice for avoiding blow-ups. 



C. Numerical examples 

In order to test the generalized harmonic shift con- 
dition we have performed a series of numerical experi- 
ments. We evolve Minkowski initial data, but with a 
non-trivial initial slice given in Minkowski coordinates 
(tM, %M) as £m = p{xm), with p a profile function that 
decays rapidly. If we use x = xm as our spatial coordi- 
nate, the spatial metric and extrinsic curvature turn out 
to be 



g(t = 0) = l-p> 2 



K(t = 0) = -p"/g . (6.25) 



In all the simulations shown below we have taken for the 
function p(x) a Gaussian centered at the origin 



pix) = k exp 



& 



(6.26) 



For our simulations we have chosen for k and s the same 
values used in |25|. namely k = 5 and s = 10. Further- 
more, we start with unit lapse and vanishing shift. 



All runs have been performed using a method of lines 
with fourth order Runge-Kutta integration in time, and 
standard second order centered differences in space. Fur- 
thermore, we have used 64,000 grid points and a grid 
spacing of Ax = 0.0125 (which places the boundaries 
at ±400), together with a time step of At = Ax/ '4. In 
the simulations shown below, we will concentrate on two 
different aspects: First, we want to know how the gen- 
eralized harmonic shift condition works in practice, and 
what are its effects on the evolution. Also, we want to 
see if gauge shocks do form when they are expected. 

Furthermore, to study the overall growth in the evo- 
lution variables, we introduce the quantity 6 defined 
through 



<5 2 := ( a -l) 2 + a 2 + (.9-l) 2 + ][>. 



2 

i ' 



(6.27) 



as a measure of how non-trivial the data is. For S we 
then also calculate the convergence factor 77 which, using 
three runs with high (S h ), medium (<5 m ) and low (5 l ) 
resolutions differing in each case by a factor of two, can 
be calculated by 



(6.28) 



In the plots we show three different convergence fac- 
tors. In particular, we denote with a triangle the conver- 
gence factor obtained when comparing runs with 64, 000, 
32, 000 and 16, 000 grid points and a spatial resolution of 
0.0125, 0.25 and 0.5. We then use boxes and diamonds 
when gradually lowering all three resolutions by a factor 
of two. For second order convergence we expect 77 ~ 4. 

As a reference of what happens for the case of zero 
shift, in Fig. [I]we show a run that corresponds to har- 
monic slicing and vanishing shift (these plots should be 
compared with Fig. 2 of j25j). In the figures, the initial 
data is shown as a dashed line, and the final values at 
t = 200 as a solid line (remember that the initial metric is 
non-trivial) . Intermediate values at intervals of At = 20 
are shown in light gray. As can be seen from the plots, 
all variables behave in a wavelike fashion and the conver- 
gence plot indicates that we have close to second order 
convergence during the whole run for all resolutions con- 
sidered. Here the pulses are moving out symmetrically 
in both directions away from the origin (we only show 
the x > side). One can see that the initial non-trivial 
distortion in g for small x remains (so the dashed and 
solid lines lie on top of each other there) , indicating that 
even though in the end we return to trivial Minkowski 
slices, we are left with non-trivial spatial coordinates. 

Our second example is shown in Fig.|21and corresponds 
to / = h = 1, that is, pure harmonic coordinates in 
both space and time. The simulation is very similar to 
the previous one and convergences again to second order. 
The non-trivial shift, however, behaves in such a way 
that at the end of the run no distortion remains in the 
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FIG. 1: For a simulation with harmonic slicing (/ = 1) and 
vanishing shift, we show the evolution in time of the variables 
a, g and K, together with that of the eigenfield urf. The 
values of the different quantities are shown every At = 20. 
In the plot on the bottom we show three convergence factors 
when increasing the resolution in the order "diamond, box 
and triangle". In all three cases the convergence factor is 
close to the expected value of 4.. 



metric component g at the origin. One should also note 
that for / = h = 1 both eigenfields propagate with the 
same speed. However, since quadratic and mixed source 
terms in the evolution equations of both w£ and w± are 
not present, simple wavelike behavior for all variables is 
again observed. 

This example allows us to understand the main effect 
that the introduction of the generalized harmonic shift 
condition has on the evolution: It drives the spatial co- 
ordinates to a situation where no final distortion in the 
metric is present. In fact, it is not difficult to under- 
stand why this is so. From equation (|4.4|) we can see 
that the sources for the evolution of the rescaled shift 
are the derivatives of the lapse, the trace of the extrinsic 
curvature and the A l mn . As the shift condition does not 
feed back into the slicing condition (apart from a triv- 
ial shifting of the time lines), the lapse and the trace of 
the extrinsic curvature behave just as before, with pulses 
that propagate away. However, the A l mn will continue to 
drive the evolution of the shift unless they become zero. 
The behavior of the shift condition is then to drive the 
system to a situation where the A l mn vanish. In the sim- 
ple 1+1 case this is equivalent to reaching a state where 
the spatial metric itself becomes trivial. 

A third example is presented in Fig.|3| which uses again 
/ = 1, but now we take h = 1 + 3er 2 . Initially, the evolu- 
tion behaves in a very similar way to the previous case. 
At later times, however, we observe that a sharp gradi- 




50 100 150 200 250 X 



T| 4.025 
4.000 
3.975 
3.950 

3 925 



50 100 150 200 250 X 



50 100 150 200 1 



FIG. 2: For a simulation with harmonic slicing and harmonic 
shift (/ = h — 1), we show the evolution in time of a, a, g 
and K, together with that of and w+. As in the previ- 
ous figure, the bottom plot shows the convergence factors for 
different resolutions. 



ent develops in the rescaled shift a, with a corresponding 
large spike in the shift eigenfield w±. Moreover, from 
the convergence plot we see that there is a clear loss of 
convergence, and as the resolution is increased, this loss 
of convergence becomes more sharply centered around 
a specific time t ~ 150, indicating that a blow-up hap- 
pens at this time. Since in this case h is a function of cr, 
the source criteria is clearly not satisfied. The fact that 
a large spike has developed in w± therefore strengthens 
the case for the source criteria being a good indicator of 
when blow-ups can be expected. 

Our final example uses / = 1 and h = 2, and is shown 
in Fig. 21 This example is interesting as the speeds for the 
lapse and shift eigenfields are different. Concentrate first 
on the evolution of the lapse, the extrinsic curvature and 
wL . Here the evolution is essentially identical to that of 
the previous examples, converging again to second order: 
a pulse travels with roughly unit speed and behind it 
everything rapidly relaxes back to trivial values. The 
eigenfield w+, on the other hand, shows a pulse traveling 
faster, with a speed ~ \/2. It also takes considerably 
longer for the region behind this pulse to relax to trivial 
values. Finally, the metric g and rescaled shift a separate 
into two pulses traveling at the two different eigenspeeds. 
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FIG. 3: For / = 1 and h = 1 + 3<t 2 , the simulation fails 
shortly before the time t = 200 due to a sharp gradient devel- 
oping in the rescaled shift a, and a corresponding large spike 
appearing in the eigenfield w±. The bottom plot shows that 
convergence starts to be lost at t ~ 150 (notice the change of 
scale as compared to previous plots), indicating that a blow- 
up has happened at around this time. This type of behavior is 
expected in this case since the source criteria is not satisfied. 



FIG. 4: For f — 1 and h — 2, the lapse and shift eigen- 
fields travel at different speeds. The lapse, extrinsic curva- 
ture and eigenfield w+ show a pulse traveling with roughly 
unit speed, while the eigenfield tu+ shows a pulse moving with 
speed ~ V2. The metric g and rescaled shift a, on the other 
hand, separate into two pulses traveling at the two different 
eigenspeeds. 



This is to be expected, as from (|6.13f) and 1)6. 14|) we see 
that metric and shift have contributions from both types 
of eigenfields. 

For the different runs we have also studied the behavior 
of the logarithm of the root mean square (rms) of 6 over 
time. Since the behavior of the evolution turns out to de- 
pend to some extend on the initial data, and in particular 
on the sign of the Gaussian in l|6.26[) . we perform runs 
for both k = 5 and n — — 5, and then take the average of 
both runs when calculating S. For the initial data we are 
using, at time t = this yields a value log(d) fts —1.583 
for both signs of k. In Fig. [3] we plot the rms of the 
quantity S for the times t = {20,40,60,80,100}, when 
using either h = 1 and varying the (constant) value of 
/ (top panel), or using / = 1 together with different 
(again constant) values of h (bottom panel). From the 
top panel we see that / = 1 is clearly preferred. In ad- 
dition we want to point out that runs with / < 0.79 and 
/ > 1.25 crashed before reaching the time t — 100. This 
behavior is expected as we know that constant values 
of / different from one produce blow-ups. In the lower 
panel we observe that for / = 1 corresponding to har- 



monic slicing, h — 1 performs best. In addition, values 
h ~ 0.5 and h ^> 1 also seem to be preferred. One should 
note that mixed terms w±w± in the evolution equations 

of both w± and w± for these choices of h play a minor 
role since localized perturbations in these eigenfields sep- 
arate quickly when traveling with different speeds. We 
also want to mention that for h < 0.19 the simulations 
again crashed before reaching the time t = 100. The ob- 
servation that S grows rapidly and runs crash early if / 
and/or h are very close to zero can be understood by the 
fact that the system is not strongly hyperbolic if / = 
and/or h = 0. 

In the contour plot of Fig. we show the rms of 6 
at time t = 100 as a function of the gauge parameters / 
and h, using 64 x 80 equidistant parameter choices. Cases 
that have already crashed by that time correspond to the 
hashed regions. Note that the darker regions in this plot 
denote parameter choices where a significant growth in 
the evolution variables is present, while brighter regions 
correspond to runs with very little growth. We find small 
values for the rms of 6 for / being close to its shock 
avoiding value / = 1, and either ft = 1 or ft > 1. In 
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FIG. 5: Top. For evolutions with h = 1, the rms of S is shown 
on a logarithmic scale as a function of / every At = 20. The 
value / = 1 is obviously preferred. Bottom. For runs with 
harmonic slicing (/ = 1), the same quantity is plotted as a 
function of h. Here h = 1 is the optimal choice, but h ~ 0.5 
or h ^> 1 is also preferred. 
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FIG. 6: Contour plot of the rms of S at time t = 100. Small 
values for this quantity are found for / = 1, when h = 1 or 
h ^> 1, and also for / = h. 



addition, we can also observe that / = h corresponds to 
a preferred choice. This can be explained by the fact that 
for this gauge choice the mixed terms w^w^. are missing 
in the evolution equations of both u>£ and w± . 



VII. EINSTEIN EQUATIONS IN SPHERICAL 
SYMMETRY 



As a second application of the generalized harmonic 
shift condition we will consider vacuum general relativ- 
ity in spherical symmetry. This situation is considerably 
richer than the 1+1 dimensional case, but it also presents 
some special problems because of the singular nature of 
spherical coordinates at the origin. 



A. ADM evolution equations 

We will consider the spherically symmetric line element 
written in the form 



ds 1 = -of i 
+Adr 2 



I- Act 2 ) dt 2 



2aAo~drdt 



Br 2 dQ 2 



(7.1) 



where all the metric coefficients are functions of both t 
and r. We now introduce the following auxiliary variables 

D a := d r hia , d a :— d r a , (7-2) 

D A :=d r ]xiA, D B :=d r \nB. (7.3) 

Notice again that we use logarithmic derivatives for the 
lapse and the spatial metric, but only an ordinary deriva- 
tive for the shift. For the extrinsic curvature, we will use 
the mixed components 



K A := K r r 



(7.4) 



Following j3l| . we will change our main evolution vari- 
ables and make use of the "anti-trace" of the metric spa- 
tial derivatives D = Da — 2Db, and the trace of the ex- 
trinsic curvature K = Ka + 2Kb, instead of Da and Ka- 
For the regularization of the evolution equations at 
the origin we will follow the procedure described in |3l| . 
which requires the introduction of an auxiliary variable 



A = (1- A/B)/r 



(7.5) 



Local flatness guarantees that A is regular and of order r 
near the origin. By taking {a, A, B,d a , K, Kb} as even 
functions at r = 0, and {u, D a , D, Db, A} as odd, one 
obtains regular evolution equations at r = 0. 

In terms of the variables introduced above, the hamil- 
tonian and momentum constraints become (in vacuum) 



= C, = -i), D L! + ^ (D + ^ 



Db 
2 



AK B (2K - 3K B ) + -{D-D B -\) 



= C m = -d r K B + {K-3K B ) 



Db 
2 



(7.6) 



(7.7) 
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Notice that the hamiltonian constraint is regular, 
while the momentum constraint still has the term 
(K — 2>K B )/r = (Ka — K B )/r which has to be handled 
with care numerically. This is not a problem as the mo- 
mentum constraint does not feed back into the ADM 
evolution equations. On the other hand, when one 
adds multiples of the momentum constraint to the evolu- 
tion equations in order to obtain strongly hyperbolic re- 
formulations (as in the following section), the regulariza- 
tion procedure requires some of the dynamical variables 
to be redefined by adding to them a term proportional 
to A (see |3lJ for details). This redefinition, however, 
does not affect the characteristic structure of the system. 
Because of this, in the following analysis we will simply 
ignore this issue. 

For the evolution of the lapse we will again take the 
Bona-Masso slicing condition, which in spherical symme- 
try takes the form 

8 t a = a 2 (aD a - fK) . (7.8) 

For the shift we will use the generalized harmonic shift 
condition in the form (|4.4|l . In this case one finds 

D 2 



7 



ry mn (3)p 

n (3)pr 



Hal. 



2A 

2 

~7b 



rA 



(7.9) 
(7.10) 



Notice that in the first of these expressions we have 
used the Christoffel symbols for the full spatial metric 
dl 2 = Adr 2 + Br 2 dfl 2 , while in the second we used those 
of the flat metric dl 2 = dr 2 + r 2 dfl 2 . However, as is clear 
from H4.4fl . in both cases we have to contract indices us- 
ing the full inverse metric which explains why there is 
a factor B in the denominator of the second expression. 
Using these expressions we then find 



Ar= D__2X 
2A A 



(7.11) 



with A defined in (|7.5|l above. Our final shift condition 
is then regular at the origin and has the form 



dt<J = a 



D a ( D 2A 
° d °- — + h {2A +aK -^ 



(7.1z 



It is important to mention that if we had used the origi- 
nal condition (|3.7I) instead of l|4.4l) , we would have found 
that the shift evolution equation was singular. More- 
over, one also finds that taking a and equal to zero 
is consistent when using (|4.4|l in the sense that their re- 
spective evolution equations guarantee that they remain 
zero, which would not have been the case with l|3.7[l . 

Going back to the metric components A and B, we 
find for their evolution equations 



d t A 



2aA 



D 
~2~ 



9 t B = 2aB 



d a -K + 2K B 
Db , 



-)-K B 

r 



(7.13) 
(7.14) 



The evolution equations for D a , d CT , D and D B again 
follow trivially from the above equations. Finally, the 
ADM evolution equations for the extrinsic curvature 
components turn out to be 



d t K = - | - d r D a - 2 d r D B + a A d r K 
+ D a (^- D a j +D b (d + ^- 



AK - - {D a - D + D B + A) 
r 



(7.15) 
D a D B , DD B 



+ AKK B - i U> a - y +D B + X ) \ . (7. 



Notice that these are directly the standard ADM evo- 
lution equations written in terms of {K, K B }, with no 
multiples of the constraints added to them. In the next 
section we will consider how such adjustments affect the 
hyperbolicity of the full system. 



B. Adjustments and hyperbolicity 

In order to analyze the characteristic structure of the 
full system of evolution equations including the gauge 
conditions, we start by defining 



u := (a,a,A,B,\) , 

v := (D a ,d a ,D,D B ,K,K B ) 



(7.17) 
(7.18) 



The system of equations can then be written in the form 
(|5.1() - (|5.2() . It turns out that by doing this, one finds 
that the ADM evolution system introduced above is not 
strongly hyperbolic when / = 1 and/or h = 1. This is 
undesirable, as these cases correspond precisely to purely 
harmonic coordinates. 

Following |26j. in order to obtain strongly hyperbolic 
systems we will consider adjustments to the evolution 
equations of the extrinsic curvature components K and 
K B of the form 



d t Vi + ^2 Aijd r Vj + hi-^ C h 

3=1 



(7.19) 



Note that we are considering only very restricted ad- 
justments here. In particular, we do not modify the evo- 
lution equations for the -D's and for d a . As explained in 
Ref. pfj . this is important for the blow-up analysis in 
the next section, as otherwise the constraints that link 
the D's to derivatives of the u's will fail to hold and the 
analysis breaks down j34[. Furthermore, for simplicity 
we will not consider adjustments that use the momen- 
tum constraint. 
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For the coefficients Hk and Kk b we make the following 
ansatz 

h K = -2 + b(a,a,A,B) , (7.20) 
h KB = [c{a,a,A,B) - l]/2 . (7.21) 

With these adjustments we find that the characteristic 
matrix for our system of evolution equations becomes 



/ -a /ON 

1/A -o- -/i/2vl -/w 

-2<7 -2 -cr 2 -8 

-a 2 

1/A b/A -a 

V c/2A -a- J 



A = 



One may now readily verify that this matrix has the fol- 
lowing eigenvalues 



a£ = o i - - 



(7.23) 
(7.24) 
(7.25) 



The system is therefore hyperbolic for {/, h, c} > 0. Fur- 
thermore, there exists a complete set of eigenvectors as 
long as c ^ / and c ^ h, so the system is strongly hyper- 
bolic except in those two cases. The eigenfields turn out 
to be: 



w±. For the quadratic source terms we find 



Jff 

-±±± 



hhh 
-±±± 



(c-f) 

1 



2 



9/1 



(7.29) 
(7.30) 



(c — /i) <9(7 

Demanding now that these terms vanish we obtain pre- 
cisely the same conditions on / and h as in the 1+1 di- 
mensional case. So again / = 1 + const/a 2 and h = h(a) 
are shock avoiding solutions. Furthermore, if one chooses 
/ = h = 1 + const/a 2 , mixed terms of the form 



(7.22) do not appear in the evolution equations of the gauge 
eigenfields w£ and w±, so this is a preferred choice. 

In contrast to the 1+1 dimensional case, now also blow- 
ups associated with the constraint eigenpair w± can arise. 
The quadratic coefficient in this case takes the form 



4±± « (i 



46 + 3c) 



(7.31) 



and by asking for this coefficient to vanish we find 

6=(l + 3c)/4. (7.32) 

From (|7.20|) and (|7.21l) we then infer that hx and Hk b 
are related by 

t -, . ?>h KB 1 
= -H x — , h KB > -- , 



(7.33) 



which is the precisely the same constraint shock avoid- 
ing half-line in the {Hk, h-K B } parameter space that was 
found in Ref. 261. 



= (c-f)D a -bfD 1 



± y/fA[(c-f)K-2bK B ] , 



(7.26) 



(c - h) A 1 ' 2 



± 



Vh &( 



1 ± aVhA K 



1 ± aVhA 



2c 



D B 



2v A 



6 1 + aVhA - 2h 



K 



B ■ 



= J~cD B ±2VAK L 



(7.27) 
(7.28) 



It is clear from these expressions that when c = / the 
first and third pairs of eigenfields become proportional to 
each other and are hence no longer independent, while 
for c = h it is the second and third pairs that become 
proportional. 



C. Gauge and constraint shocks 

As we did for the 1+1 dimensional system, we will now 
study the possible formation of blow-ups for the evolu- 
tion equations in spherical symmetry. In order to apply 
the source criteria for avoiding blow-ups we need to calcu- 
late the quadratic source terms in the evolution equations 
for the eigenfields. We first look for "gauge shocks" , for 
which we concentrate on the gauge eigenfields w± and 



D. Numerical examples 

We will now test the effects of the generalized harmonic 
shift in spherical symmetry by performing a series of nu- 
merical simulations. As in the 1+1 dimensional case, we 
will concentrate on two aspects, namely the effect of the 
shift condition on the evolution of the metric, and the 
possible formation of blow-ups. Furthermore, in order to 
decouple geometric effects associated with the center of 
symmetry, we will consider two distinct regimes, one far 
from the origin and one close to it. 



1. Pulses far from the origin 

We will first consider simulations that are far from the 
origin, using initial data that is similar to the one used in 
Sec. IVI CI We start with the Minkowski spacetime, but 
use a non-trivial initial slice with a profile <m = p(tm)- 
The initial metric and extrinsic curvature then become 



A(t = 0) 
B(t = 0) 

K(t = 0) 
K B (t = 0) 



1-P'\ 
1 , 

VA\A r 

P' 
rs/A 



(7.34) 
(7.35) 

(7.36) 
(7.37) 
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The profile function p(r) is again chosen to be a Gaussian, 



p(r) — k exp 



(7.38) 



using for its amplitude and width the values k = ±5 and 
s = 10. The center of the Gaussian is taken at r c = 250, 
such that for evolution times of t ~ 100 the perturbation 
will remain away from the origin. 

We have performed runs with the code described 
in jsJl, which uses a method of lines with fourth order 
Runge-Kutta integration in time, and standard second 
order centered differences in space. We used 5,000 grid 
points and a grid spacing of Ar — 0.1 (which places 
the outer boundary at 500) together with a time step 
of At = Ar/4. 

In a preparatory experiment, we studied which evolu- 
tion systems perform best for our optimal gauge choice 
f = h = 1 . The upper panel of Fig. (which should be 
compared with Fig. 7 of 26] ) shows the rms of the hamil- 
tonian constraint at time t = 100, as a function of the 
adjustment parameters hx and Hk b ■ We can see that the 
line H7.33f) obtained by the source criteria is numerically 
preferred, although there does seem to be a discrepancy 
for large values of Hk b f° r which the numerical results 
suggest a somewhat steeper line. This discrepancy is due 
to the effect of 1/r terms which are not taken into ac- 
count by the source criteria and can be eliminated by 
removing these terms by hand. It is also important to 
point out that, in contrast to Ref. [2||, the initial data 
used here satisfies the constraints and all subsequent con- 
straint violations are caused by truncation error. 

In order to determine which points on this line perform 



best, 



to fix the eigenspeeds A5_ of the constraint 



mode, we tested different (constant) values of c. From the 
lower panel of Fig. [7] (to be compared with Fig. 5 of |26j) 
we find that values c ~ 1/4 and c> 1 are preferred. This 
observation can be readily understood by the fact that 
the system is not strongly hyperbolic for c = and c = 1 , 
and by the fact that for c ~ 1 we expect contributions 

and w%. 



from mixed source terms, since then w 
propagate with similar or even identical eigenspeeds. 

For our main experiment regarding gauge effects, we 
concentrated on evolution systems which belong to the 
shock avoiding family = — 1 + 3/iif B /2, where for c 
we considered three different values: c = {1/4, 1,4}. As 
long as the pulses remain far from the origin, we have 
found that the evolutions behave in a very similar way to 
those of the 1+1 dimensional case described in Sec. I VI CI 
We summarize these results in Fig. [8] showing for these 
three choices of c the rms of the hamiltonian constraint at 
time t — 100 as a function of / and h. These graphs are 
very similar to Fig. and show that / = 1 together with 
h = 1 or h 3> 1 , and f = h are again preferred parameter 
choices, indicating that the same mechanisms as in the 
1+1 dimensional case are at work. One should observe 
the different scales when comparing the three plots corre- 
sponding to different values of c, which indicate that by 
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FIG. 7: Top. Contour plot of the rms of the hamiltonian 
constraint at time t — 100 as a function of the adjustment 
parameters Hk and hx B ■ The parameter line suggested by 
the source criteria, hx = —1 + 3fiK B /2 with Yik b > —1/2, is 
shown as a solid line. Bottom. For reasons described in the 
text, along this line c ~ 1/4 and c> 1 are preferred values for 
c, the latter determining the eigenspeeds AJ. of the constraint 
modes. 



far the lowest constraint violations are found when the 
constraint eigenspeed is different from the gauge eigen- 
speeds. Notice also in the middle plot corresponding to 
c = 1 that the region around / = h = 1 is in fact dark. 
This can be explained by the fact that for / = h = c = 1 
the evolution system is not strongly hyperbolic. 

When the pulses come close to the origin, however, 
additional effects arise due to 1/r terms. In the next 
section we will consider this situation. 



2. Pulses close to the origin 

In order to see directly the effect of the generalized har- 
monic shift condition on the evolution of the geometric 
variables, we will consider again a series of simulations 
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FIG. 8: Contour plot of the rms of the hamiltonian constraint at time t — 100 for c = {1/4, 1,4}. As in the 1+1 dimensional 
case, / = 1 together with h = 1 or h 3> 1, and f — h are preferred parameter choices. Furthermore, by far the lowest constraint 
violations are found when the constraint eigenspeed is different from the gauge eigenspeeds, i.e. for 1. 



of Minkowski spacetime, but this time close to the origin 
r = 0. The initial data for these runs will be simpler than 
the one used in the previous section: We start with a flat 
Minkowski slice with A = B = 1 and Ka = Kb = 0, and 
take a non-trivial initial lapse of the form 

a = l + Kr 2 (e-^-^ 2 l s2 +e-^+ r ^' s2 ) , (7.39) 

with k = 10 -5 , r c = 10 and s — 1 (the reason for the 
two gaussians is to make sure the initial lapse is an even 
function of r). All simulations shown here use 4,000 grid 
points, with a grid spacing of Ar = 0.01 (which places 
the outer boundary at r = 40), together with a time 
step of Ar/4. In the plots, the initial data is shown as 
a dashed line and the final values at t = 20 as a solid 
line. Intermediate values are plotted every At = 2 in 
light gray. 

As reference, we first show in Fig. [5] a run for the case 
of harmonic slicing (/ = 1) with no shift. In order to 
look at the details in a clearer way, in the figure we plot 
a — 1, A — 1 and B — 1. As expected, the perturbation 
pulse in the lapse separates into two pulses, one moving 
outward and one inward. The inward moving pulse goes 
through the origin and starts moving out much in the 
way a simple scalar wave would. The pulses in the lapse 
are accompanied by similar pulses in the metric variables 
A and B. However, one can clearly see that the metric 
variables are not evolving toward trivial values, so in the 
end we are left with Minkowski slices with non-trivial 
spatial coordinates. 

Next we consider the same situation, but now using a 
harmonic shift with h = 1. Fig.llOlshows results from this 
run. The lapse behaves in exactly the same way as before, 
but now there is a non-trivial shift. The evolution of a 
indicates that the shift behaves much in the same way as 
the lapse, with two pulses traveling in opposite directions, 
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FIG. 9: Evolution of pulses close to the origin for harmonic 
slicing (/ = 1) and zero shift. In order to make the details 
more visible, here we plot a — 1, A — 1 and B — 1. 



with the inward moving pulse going through the origin 
and then moving out as expected. The evolution of the 
metric variables A and B shows that after the ingoing 
pulse goes through the origin and starts moving out, the 
perturbations on the metric become very small. The shift 
then seems to be having a similar effect to the one it had 
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FIG. 10: As previous plot, but using the harmonic shift con- 
dition with h — 1. 



in the 1+1 case, making the metric components evolve 
toward trivial values. 

Fig. II II shows a similar run, but now using / = 1 and 
h = 2. The whole simulation behaves much the same way 
as before, except for the fact that the metric coefficient 
A (and to a lesser extent B) now shows evidence of two 
pulses separating and traveling at different speeds after 
the rebound through the origin. 

Finally, in Fig. ^] we show a simulation with h — 1 for 
a case where we have left the lapse equal to one through- 
out the evolution. The initial data in this case is purely 
Minkowski data with a shift of the form 



2 /s 2 



-(r+r c ) 2 /s ; 



(7.40) 



with k = 10~ 3 , s = 1 and r c = 5. The purpose of this run 
is to decouple the harmonic shift condition from the slic- 
ing condition. The figure shows clearly how, even though 
the initial pulse in the shift produces perturbations in the 
metric coefficients, these perturbations rapidly decrease 
in size leaving trivial values behind them. 

One should mention the fact that, even though we 



FIG. 11: Same as previous plots, but now using the general- 
ized harmonic shift with h = 2. The metric coefficient A (and 
to a lesser extent B) shows evidence of two pulses separating 
and traveling at different speeds. 



don't show convergence plots in this section, in all cases 
convergence has been studied and we have found that the 
simulations converge at close to second order. 



VIII. DISCUSSION 

We have proposed a natural generalization of the con- 
dition for harmonic spatial coordinates analogous to the 
generalization of harmonic time slices of Bona et al. |19| , 
and closely related to shift conditions recently introduced 
by Lindblom and Scheel [g, and by Bona and Palen- 
zuela p^ . This coordinate condition implies an evolu- 
tion equation for the shift components. We have also 
found that if one wants to decouple this evolution equa- 
tion for the shift from the choice of slicing condition, it is 
important to work with a rescaled shift vector a 1 — f3 l /a. 

The generalized harmonic shift condition thus obtained 
turns out not to be 3-covariant, which is not surprising 
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FIG. 12: Trivial Minkowski slices where the lapse remains 
equal to unity, together with the harmonic shift condition 
(h = l). 



as it involves the 3-Christoffel symbols directly. In order 
to be able to use this condition in arbitrary sets of curvi- 
linear coordinates, and to be sure that we always obtain 
the same shift independently of the choice of spatial co- 
ordinates, we have proposed that the condition should 
be interpreted as always being applied to topologically 
Cartesian coordinates, and later rewritten in a general 
curvilinear coordinate system. In this way we have ob- 
tained a fully 3-covariant version of the generalized har- 
monic shift condition. 

We have shown that the evolution equation for the shift 
proposed here can be seen to lead to strongly hyperbolic 
evolution systems both in the case of 1+1 "toy" relativ- 
ity and in the case of spherical symmetry. Though we 
have not done a completely general analysis here, it is 
to be expected that it will also lead to strongly hyper- 
bolic systems in the 3D case. Here we have concentrated 
on simple one-dimensional systems in order to take the 
hyperbolicity analysis further and study the possible for- 
mation of blow-ups associated with this shift condition. 
We find that the coefficient h controlling the gauge speed 
associated with the shift can be an arbitrary function of 
the lapse, but must be independent of the shift itself in 
order to avoid blow-ups. In the slicing and constraint 
sectors we recover previous results found in |26| . An im- 
portant result of this study is the fact that evolutions 
will be much better behaved if the gauge speeds associ- 
ated with the lapse and shift are the same. This can be 
understood from the fact that terms in the sources that 



are mixed products of eigenfields associated to the lapse 
and shift vanish in this case. This implies that if one 
wants to use a shift of the generalized harmonic family, 
together with a lapse of the Bona-Masso type (like 1+log 
slicing), it is best to take h{a) = f{ot). In particular, the 
shock avoiding family h = f = 1 + const/a 2 is an optimal 
choice. 

We have also performed a series of numerical simu- 
lation both to confirm the predictions of the blow-up 
analysis, and to study what effect the shift has on the 
evolution of the geometric variables. In the 1+1 dimen- 
sional case, we find that the effect of the shift is to take 
the spatial metric back to a trivial value everywhere, by 
propagating away any non-trivial values in a wavelike 
fashion. In spherical symmetry the situation is consid- 
erably richer, but our main result is that when one uses 
the 3-covariant version of the generalized harmonic shift, 
then the effect of the shift is also to drive the metric coef- 
ficients to trivial values by propagating away any initial 
perturbations in the way one would expect for spheri- 
cal waves, i.e. the perturbations become smaller as they 
propagate outwards. It is important to mention that, had 
we not used the covariant form of the shift condition and 
tried to apply the original non-covariant version directly 
to spherical coordinates, we would have found the shift 
condition to be singular, and worse still, to break the 
original spherical symmetry of the system. This shows 
that working with the 3-covariant version is the correct 
approach. 

As a final comment, one should also mention the fact 
that the requirement of 3-covariance is not satisfied by 
some recently proposed shift conditions that are currently 
being used by large scale 3D simulations, such as the 
"Gamma driver" shift 0, 13^, HI] . We are currently also 
studying 3-covariant versions of those conditions. 
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APPENDIX A: GENERALIZED HARMONIC 
LAPSE AND SHIFT CONDITIONS 

Here we will provide a general derivation of equa- 
tions H3.1f) and (|3.6|l for the lapse and shift. Let us start 
by considering the dAlambertian of any number a of 
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functions ip a (x^) with their corresponding source terms 
n^ a = S r . (Al) 
Now, the d'Alambertian can be written in general as 
1 



□v>° 



-9 



(A2) 



Using g^ v — 7 A " / — n^n" , with 7 A " / the projector operator 
on the hypersurfaces S t with normal n M , we find 



1 

a^7 



d^a^n^d^} , (A3) 



where we used the fact that g := det g^ v — —a 2 j with 
7 := det 7ij the determinant of the 3-metric on Sj. We 
then have 



n"V M (r^V^ a ) , (A4) 



where 3 A is the Laplacian compatible with the 3-metric 
a M = n "V„n" = 7^V„[lna] =: D^[ln a] is the 
4-acceleration of the normal observers, and we used 
K = -V v n v . 

In order to obtain for instance a system of first order 
equations one can further define 



D"iP a , (A5) 
C n V> a = n v V y ^ a , (A6) 



where D^if)" 1 :— 7 A " T V <T i/' a - Collecting the above results 
we obtain 



C H IF - a^ a - D^ a - U a K = -S^' 



(A7) 



A simple application of the above results is the case 
when ip a — x a , in which case IF = n a = (1,— (3 l )/ct, 
Ql = 7* and D u Q vi =dj 0/77*) A/7 = - (3) ^- The 
above equation with = is then called the harmonic 
coordinate condition, which provides an evolution equa- 
tion for the lapse - Eq. (|2.9(1 - and for the shift - Eq. (|2.10|) 
- when taking x a — (t, x l ) with t defining the time slicings 
and x % being spatial coordinates on E f . 

On the other hand one can take dijj a = with a 
source term of the form = q^an^n^V M V !/ -i/ ,a ( n0 sum 
over index a). Now, using n^n u — 7^" — g^ , one obtains 



(A8) 



Using the orthogonal decomposition V ' u ij} a = D u if: a — 
n„n a V a ?P a =Ql- n u Il a we find 

n'WV M V \i\) a = D v Q va + Ii a K - Oip a , (A9) 



where we used 7 At<T V M n CT = V M n M = —K and 7 Miy n„ = 0. 
In this way the equation Dtp a = becomes 



1 + Qty< 



(B v Q va + U a K) 



(A10) 



Finally, —□■0° is given by the left hand side of 
Eq. i|A7(l . from where we find 

C H W - a^ a - D^ a - IPX = 



1 + g^< 



which simplifies to 



(A,Q" a + n a if) , (Ail) 



(D v Q va + U a K) . (A12) 



1 + qty" 

In this way by taking tjj a — (t 7 x l ), q t =at = 1/f — 1, 
q x i = ah = l/h—1, together with Eqs. I|A5[) and (|A6() 
(leading to IP = n a = (1,-0*) /a and Qg = 7°), one 
recovers the evolution equations l|3.1ll and (|3.ti[) for a and 
respectively. 



APPENDIX B: METRIC AND CHRISTOFFEL 
SYMBOLS IN 3+1 LANGUAGE 

For the expression of the generalized harmonic gauge 
conditions one needs to write the 4-metric of spacetime 
and its associated Christoffel symbols in 3+1 language. 
The 4-metric in terms of 3+1 quantities has the form 



5oo = ~(a 2 -HjPP) 
goi = Hj P J , 

Qij = lij , 

and the corresponding inverse metric is 
,00 



<T = -l/a 
9 



9 



(Bl) 
(B2) 
(B3) 



(B4) 
(B5) 
(B6) 



From this one can obtain the following expressions for 
the 4-Christoffel symbols in terms of 3+1 quantities 

r° n = (d t a + /3 m d m a - f3 m p n K mn ) /a , (B7) 

(B8) 
(B9) 



00 
^0 



r& - (dta - f3 m K m ) /a , 

r°- - -Kij/a, 

T l 00 = ad'a-2a/3 m K l m 

- P l (d t a + P m 3 m a - (3 m (3 n K mn ) /a 

+ d t f3 l +l3 m (3) V m /3 z , 

rLo = -p l (d m a-f3 n K mn )/a 

- aK l m + ^V m /3 l , 



(BIO) 

(BH) 
(B12) 
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The contracted Christoffel symbols T a := g^T% v then and 
become 

r = ^(d t a- l3 m d m a + a 2 K) + ^r 

r° = + (bu) _ | (a/3 ._ wl + aa . a) . (B14) 
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